arXiv:1407.1686v2 [astro-ph.HE] 6 Oct 2014 


Milagro Limits and HAWC Sensitivity for the 
Rate-Density of Evaporating Primordial Black Holes 


A. A. Abdo a ’ a t A. U. Abeysekara a , R. Alfaro b , B. T. Allen af ’ ak , C. 
Alvarez c , J. D. Alvarez d , R. Arceo c , J. C. Arteaga-Velazquez d , T. 
Aune aa,al , H. A. Ayala Solares 6 , A. S. Barbed, B. M. Baughman®, N. 
Bautista-Elivar h , J. Becerra Gonzalez 1 ’®, E. Belmont b , S. Y. BenZvi ae , D. 

Berley®, M. Bonilla Rosales k , J. Braun®, R. A. Caballero-Lopez 1 , K. S. 
Caballero-Mora m , A. Carraminana k , M. Castillo 11 , G. E. Christopher a ®, U. 
Cotti d , J. Cotzomi 11 , E. de la Fuente 0 , C. De Leon d , T. DeYoung a,p , R. 
Diaz Hernandez k , L. Diaz-Cruz 11 , J. C. Diaz-Velez-i, B. L. Dingus q , M. A. 
DuVernois' , R. W. Ellsworth 1 ’®, D. W. Fiorinh , N. Fraija s , A. Galindo k , F. 
Garfias s , M. M. Gonzalez s , J. A. Goodman®, V. Grabski b , M. Gussert 1 , Z. 
Hampel- Arias-i , J. P. Harding q , E. Hays 1 , C. M. Hoffman q , C. M. Hui e , P. 
Hiintemeyer 6 , A. Imran' . A. Iriarte s , P. Kara' , D. Kieda f , B. E. 
Kolterman ag , G. J. Kunde q , A. Lara 1 , R. J. Lauer 11 , W. H. Lee s , D. 
Lennarz v , H. Leon Vargas b , E. C. Linares d , J. T. Linnemann a , M. Longo 1 , 
R. Luna-Garcla w , J. H. MacGibbon ad , A. Marinelli b , S. S. Marinelli a , H. 
Martinez™, O. Martinez 11 , J. Martmez-Castro w , J. A. J. Matthews 11 , J. 
McEnery 1 , E. Mendoza Torres k , A. I. Mincer a ®, P. Miranda- Romagnoli x , E. 

Moreno 11 , T. Morgan ah , M. Mostafa p , L. Nellen y , P. Nemethy 3 ®, M. 
Newbold f , R. Noriega-Papaqui x , T. Oceguera-Becerra°’ b , B. Patricelli s , R. 
Pelayo w , E. G. Perez-Perez h , J. Pretz p , C. Riviere g,s , D. Rosa-Gonzalez k , 
E. Ruiz-Velasco b , J. Ryan 2 , H. Salazar 11 , F. Salesa p , A. Sandoval 11 , P. M. 
Saz Parkinson aa,ai , M. Schneider aa , S. Silich k , G. Sinnis q , A. J. Smith®, D. 
Stump a , K. Sparks Woodle p , R. W. Springer 1 , I. Taboada v , P. A. Toale ab , 
K. Tollefson a , I. Torres k , T. N. Ukwatta a,q , V. Vasileiou®’ am , L. Villasenor d , 
T. Weisgarbed, S. WesterhofP, D. A. Williams 33 , I. G. Wished, J. Wood®, 
G. B. Yodh ac , P. W. Younk q , D. Zaborov p , A. Zepeda™, H. Zhou 6 

° Department of Physics and Astronomy, Michigan State University, East Lansing, MI, 

USA 

b Instituto de Fisica, Universidad Nacional Autonoma de Mexico, Mexico D.F., Mexico 
c CEFyMAP, Universidad Autonoma de Chiapas, Tuxtla Gutierrez, Chiapas, Mexico 
d Universidad Michoacana de San Nicolas de Hidalgo, Morelia, Mexico 
e Department of Physics, Michigan Technological University, Houghton, MI, USA 
1 Department of Physics and Astronomy, University of Utah, Salt Lake City, UT, USA 
9 Department of Physics, University of Maryland, College Park, MD, USA 
h Universidad Politecnica de Pachuca, Pachuca, Hgo, Mexico 
'NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA 
■'Department of Physics, University of Wisconsin-Madison, Madison, WI, USA 


1 



k Instituto Nacional de Astrofisica, Optica y Electronica, Tonantzintla, Puebla, Mexico 
l Instituto de Geofisica, Universidad Nacional Autonoma de Mexico, Mexico D.F., Mexico 
m Physics Department, Centro de Investigacion y de Estudios Avanzados del IPN, Mexico 

City, DF, Mexico 

n Facultad de Ciencias Fisico Matematicas, Benemerita Universidad Autonoma de 

Puebla, Puebla, Mexico 

° IAM-Dpto. de Fisica; Dpto. de Electronica (CUCEI), IT.Phd (CUCEA), Phys-Mat. 

Phd (CUVALLES), Universidad de Guadalajara, Jalisco, Mexico 
p Department of Physics, Pennsylvania State University, University Park, PA, USA 
q Physics Division, Los Alamos National Laboratory, Los Alamos, NM, USA 
r School of Physics, Astronomy, and Computational Sciences, George Mason University, 

Fairfax, VA, USA 

s Instituto de Astronomia, Universidad Nacional Autonoma de Mexico, Mexico D.F., 

Mexico 

t Colorado State University, Physics Dept., Ft Collins, CO 80523, USA 
u Dept of Physics and Astronomy, University of New Mexico, Albuquerque, NM, USA 
v School of Physics and Center for Relativistic Astrophysics - Georgia Institute of 
Technology, Atlanta, GA, USA 30332 

w Centro de Investigacion en Computacion, Instituto Politecnico Nacional, Mexico City, 

Mexico 

x Universidad Autonoma del Estado de Hidalgo, Pachuca, Hidalgo, Mexico 
v Instituto de Ciencias Nucleares, Universidad Nacional Autonoma de Mexico, Mexico 

D.F., Mexico 

z Space Science Center, University of New Hampshire, Durham, NH, USA 
aa Santa Cruz Institute for Particle Physics, University of California, Santa Cruz, Santa 

Cruz, CA, USA 

ab Department of Physics & Astronomy, University of Alabama, Tuscaloosa, AL, USA 
ac Department of Physics and Astronomy, University of California, Irvine, Irvine, CA, 

USA 

ad Department of Physics, University of North Florida, Jacksonville, FL 32224, USA. 
ae Department of Physics and Astronomy, University of Rochester, Rochester, NY, USA. 
a f Department of Physics and Astronomy, University of California, Irvine, CA 92697 
Department of Physics, New York University, 4 Washington Place, New York, NY 

10003 

ah Department of Physics, University of New Hampshire, Morse Hall, Durham, NH 03824 
“ Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong, 

China 

a] Current address: Operational Evaluation Division, Institute for Defense Analyses, 4850 
Mark Center Drive, Alexandria, VA 22311-1882 
ak Current address: Harvard- Smithsonian Center for Astrophysics, Cambridge, MA 02138 
al Current address: Department of Physics and Astronomy, University of California, Los 

Angeles, CA 90095 

am Current address: Laboratoire Univers et Particules de Montpellier, Universite 
Montpellier 2, CNRS/IN2P3, CC 72, Place Eugene Bataillon, F-34095 Montpellier 

Cedex 5, France 


2 



Abstract 


Primordial Black Holes (PBHs) are gravitationally collapsed objects that 
may have been created by density fluctuations in the early universe and 
could have arbitrarily small masses down to the Planck scale. Hawking 
showed that due to quantum effects, a black hole has a temperature in- 
versely proportional to its mass and will emit all species of fundamental 
particles thermally. PBHs with initial masses of ~ 5.0 x 10 14 g should be 
expiring in the present epoch with bursts of high-energy particles, including 
gamma radiation in the GeV - TeV energy range. The Milagro high energy 
observatory, which operated from 2000 to 2008, is sensitive to the high end 
of the PBH evaporation gamma-ray spectrum. Due to its large field-of-view, 
more than 90% duty cycle and sensitivity up to 100 TeV gamma rays, the 
Milagro observatory is well suited to perform a search for PBH bursts. Based 
on a search on the Milagro data, we report new PBH burst rate density up- 
per limits over a range of PBH observation times. In addition, we report 
the sensitivity of the Milagro successor, the High Altitude Water Cherenkov 
(HAWC) observatory, to PBH evaporation events. 

Keywords: Primordial Black Holes 


1. Introduction 

Primordial Black Holes (PBHs) are created from density inhomogeneities 
in many scenarios of the early universe [lj. The initial mass of a PBH is 
expected to be roughly equal to or smaller than the horizon or Hubble mass 
at formation, giving possible PBH masses ranging from that of supermas- 
sive black holes down to the Planck mass. PBH production can thus have 
observable consequences today spanning from the largest scales, for example 
influencing the development of large-scale structure in the Universe, to the 
smallest scales, for example enhancing local dark matter clustering. Addi- 
tionally, PBHs in certain mass ranges may constitute a fraction of the dark 
matter in the universe [1]. For particle physics, the greatest interest is in the 
radiation directly emitted by a black hole. By evolving an ingoing solution 
past a gravitationally collapsing object, Hawking showed that a black hole 
will thermally emit (‘evaporate’) all available species of fundamental parti- 
cles [2] with a temperature inversely proportional to the black hole mass. 
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PBHs with an initial mass of ~ 5.0 x 10 14 g should be expiring now with 
bursts of high-energy particles, including gamma radiation in the GeV - 
TeV energy range [3] . 

Detection of radiation from a PBH burst would provide valuable in- 
sights into the early universe and many areas of physics, as well as confirm 
the amalgamation of classical thermodynamics with general relativity. Ob- 
served PBH radiation will give access to particle physics at energies higher 
than those which will likely ever be achievable in terrestrial accelerators. 
Non-detection of PBHs in dedicated searches also gives important informa- 
tion. One of the most important cosmological motivations for PBH searches 
is to place limits on the initial density fluctuation spectrum of the early 
universe. In particular, PBHs can form from the quantum fluctuations asso- 
ciated with many types of inflationary scenarios |3|. Other PBH formation 
mechanisms include those associated with cosmological phase transitions, 
topological defects or an epoch of low pressure (soft equation of state) in 
the early universe (for a review see [1]). 

Evaporating PBHs are candidate sources for gamma-ray telescopes, as 
they produce short bursts of gamma rays. While most gamma-ray bursts 
(GRBs) are generally thought to be produced by the collapse of massive 
stars (long duration GRBs) or the merger of compact objects (short dura- 
tion GRBs) at cosmological distances [5], some short duration GRBs show 
behavior, such as large offsets from the host galaxy or anisotropic sky distri- 
bution, that may indicate a more local origin. Models of PBH evaporation 
based on Standard Model physics predict a unique TeV gamma-ray spec- 
trum Tjj. 

Various detectors have searched for PBH events using direct and in- 
direct methods. These methods probe the PBH distribution on various 
distance scales. One can probe the PBH density on the cosmological scale 
using the 100 MeV extragalactic ganuna-ray background, which produces 
a limit on the corresponding cosmological average PBH burst rate density 
of < 10 -6 pc~ 3 yr~ - 1 EE]. On the galactic scale, if PBHs are clustered 
in the Galaxy, we would expect to see an enhancement in the local PBH 
density and anisotropy in the 100 MeV gamma-ray measurements. Indeed, 
such an anisotropy has been detected and results in a corresponding Galac- 
tic PBH burst limit of < 0.42 pc _3 yr _1 jS]. On the kiloparsec scale, the 
Galactic antiproton background can be used to give a PBH burst limit of < 
0.0012 pc -3 yr~ 1 m- However, the antiproton-derived limit depends on the 
assumed PBH distribution within the Galaxy and the propagation of an- 
tiprotons through the Galaxy, as well as the production and propagation of 
the secondary antiproton component produced by interactions of cosmic-ray 
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nuclei with the interstellar gas. On the parsec scale, the PBH burst limits 
are directly set by searches for the detection of individual bursting PBHs 
and are independent of assumptions of PBH clustering. The best direct 
search limits come from the Very High Energy (VHE) searches conducted 
with Imaging Air Cherenkov Telescopes (IACTs) and Extensive Air Shower 
(EAS) arrays. On the parsec scale, the current best PBH limit from direct 
searches is < 1.4 x 10 4 pc _3 yr _1 [TO] , 

2. Milagro and HAWC Observatories 

In this paper, we present new PBH burst limits based on the direct 
search method using the data from the Milagro observatory. These limits 
are obtained assuming the standard model of Hawking radiation and particle 
physics eh El- Milagro was a water Cherenkov gamma-ray observatory 
(EAS type) sensitive to gamma rays in the energy range ~ 50 GeV to 100 
TeV. The observatory was located near Los Alamos, New Mexico, USA 
at latitude 35.9° north, longitude 106.7° west and an altitude of 2630 nr, 
and was operational from 2000 to 2008 |12l . The Milagro detector had two 
components: a central rectangular 60 m X 80 nr X 7 nr reservoir filled with 
purified water and an array of 175 smaller outrigger (OR) tanks distributed 
over an area of 200 nr X 200 m surrounding the reservoir. The reservoir was 
light-tight and instrumented with two layers of 8-inch photomultiplier tubes 
(PMTs). The top layer consisted of 450 PMTs (the ‘air-shower’ layer) 1.5 
nr below the water surface and the bottom layer had 273 PMTs (the ‘muon’ 
layer) 6 m below the surface. Each outrigger tank contained one PMT. The 
observatory detected VHE gamma rays by detecting the Cherenkov light 
generated by the secondary particles from the gamnra-ray-induced air shower 
as the secondary particles passed through the water. Various components 
of the detector were used to measure the direction of the gamma ray and to 
reduce the background due to hadron-induced showers. Because of its large 
field-of-view, ~2 sr, and high duty cycle, over 90%, Milagro was well suited 
to search for burst emission from PBHs. 

In this paper we also present the sensitivity of the High Altitude Water 
Cherenkov (HAWC) observatory to PBH bursts. HAWC, the successor to 
Milagro, is the next generation VHE observatory presently under construc- 
tion at Sierra Negra, Mexico at an altitude of 4100m. HAWC will consist 
of 300 water tanks, each 7.3 m in diameter and 4.5 m deep. Each tank will 
house three 8-inch PMTs (reused from Milagro) and one 10-inch PMT [13] . 
The PMTs will detect Cherenkov light from secondary particles created in 
extensive air showers induced by VHE gamma rays of energy in the range 
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~50 GeV to 100 TeV. HAWC has two data acquisition (DAQ) systems: the 
main DAQ and the scaler DAQ. The main DAQ system measures the arrival 
direction and energy of the high-energy gamma-rays by timing the arrival 
of particles on the ground. The direction of the original primary particle 
may be resolvable with an error of between 0.1 and 2.0 degrees depending 
on its energy and location in the sky. The scaler DAQ system counts the 
number of hits in each PMT, allowing a search for excesses over background 
noise. The scaler DAQ system is more sensitive to lower energy air showers 
than the main DAQ system. HAWC has a large field-of-view (1.8 sr corre- 
sponding to 1/7 th of the sky) and will have a high duty cycle of greater 
than 90%. Thus, HAWC will be able to observe high-energy emission from 
gamma-ray transients [14j . 

3. Methodology 

3.1. Primordial Black Hole Burst Spectrum 

The properties of the final burst of radiation from a PBH depend on 
the physics governing the production and decay of high-energy particles. 
As the black hole evaporates, it loses mass and hence its temperature and 
the number of particle species that it emits increase until the end of its 
evaporation lifetime. In the Standard Evaporation Model (SEM) [11* H5] . 
a PBH should directly radiate those fundamental particles whose Compton 
wavelengths are on the order of the size of the black hole. When the black 
hole temperature exceeds the Quantum Chromodynamics (QCD) confine- 
ment scale (250-300 MeV), quarks and gluons should be directly emitted 
by the black hole mm- The quarks and gluons will then fragment and 
hadronize as they stream away from the black hole, analogous to the jets 
seen in terrestrial accelerators mm- On astrophysical timescales, the jet 
products will decay into photons, neutrinos, electrons, positrons, protons 
and anti-protons. 

Detailed studies using the SEM to simulate the particle spectra from 
black holes with temperatures of 1 — 100 GeV have shown that the gamma- 
ray spectrum is dominated by the photons produced by neutral pion decay 
in the Hawking-radiated QCD jets and is broadly peaked at photon energies 
of ~100 MeV. The photons which are directly radiated (not the result of 
decay of other primary particles) are visible as a much smaller peak at a 
much higher photon energy proportional to the black hole temperature HU. 
As the evaporation proceeds to higher temperatures, the greater will be 
the number of emitted fundamental particle degrees of freedom. As new 
particle degrees of freedom become available, the luminosity of the burst 
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will increase. The energy spectrum is determined by the correct high energy 
particle physics model. In this work, we will assume the SEM as our emission 
and particle physics model. 

For temperatures well below the Planck temperature, the temperature 
(T) of a black hole can be expressed in terms of the remaining evaporation 
lifetime (r) of the black hole (that is, the time left until the black hole stops 
evaporating) as follows [15) ! 16J : 


r~ 


4.8 x 10 11 



-1 1/3 

GeV. 


( 1 ) 


The emission rate increases as the black hole shrinks and therefore the re- 
maining evaporation lifetime decreases as T increases. For black holes with 
temperatures greater than several GeV at the start of the observation, the 
time integrated photon flux can be parameterized as m 


—— k 9 x 10 35 particles m 2 GeV 1 
dE 


( lGeV yl* ( lGeV y/* 

(±^) 3 , E >T 


E <T 


(2) 


for gamma ray energies E > 10 GeV. The E~ 3 fall off at E > T comes 
from the r oc T~ 3 dependence in Equation [TJ which is less steep than the 
high energy exponential tail in the instantaneous Hawking spectrum at each 
T [T6]. Figure [I] shows the resulting gamma-ray spectrum for various PBH 
remaining lifetimes ranging from 0.001 s to 100 s. 


3.2. Detectable Volume Estimation 

In order to calculate the upper limits on the PBH burst rate density, it 
is necessary to calculate the PBH detectable volume for a given detector. In 
general, the expected number of photons detectable by an observatory on 
the Earth’s surface from a PBH burst of duration r at a non-cosmological 
distance r and zenith angle 6 is 


li(r,9,r) 


(1 -/) [ E2 

4vrr 2 J El 


dN(r) 

dE 


A(E, 6) dE 


( 3 ) 


where / is the dead time fraction of the detector and dN/dE is the black 
hole gamma ray spectrum integrated from remaining time r to 0. The values 
E\ and correspond to the lower and upper bounds respectively of the 
energy range searched and A(E, 9) is the effective area of the detector as a 
function of photon energy and zenith angle. Typically the function A(E, 9) 
is obtained from a simulation of the detector. For Milagro and HAWC, 
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Figure 1: Time-integrated gamma-ray spectrum over various PBH remaining lifetimes 
using the parametrization of Equation [2] The black hole temperature at the start of 
observation is shown in parentheses. 


the dependence of A(E, 9) on the zenith angle is usually given in discrete 
bands (represented by 6j). We will define these bands specifically for a given 
observatory in Section [4j 

The background cosmic-ray flux at the Earth’s surface is much higher 
than the background gamma-ray flux. Most events detected by EAS arrays 
such as Milagro or HAWC are air showers induced by cosmic rays. To 
search for the emission from a PBH burst one needs to look for an excess 
that cannot be explained by statistical fluctuations of the background. 

In this paper, we estimate /x o (0j,r), the minimum number of counts 
needed for a detection for different burst durations r, by finding the number 
of counts required over the background for a 5 a significance (after correct- 
ing for multiple trials). First, we calculate the background rates ( R(9i )) 
utilizing a Monte Carlo simulation (for HAWC) or actual data (for Mila- 
gro). The background rate depends on the spatial bin-size, burst duration 
and background rejection parameters, as well as the zenith angle band 9^. In 
section [4j we optimize these parameters to minimize the background rates 
and maximize the sensitivity. Using the background rates, we then find the 
Ho (9i, t) values required for a 50% probability of detecting a 5cr excess after 
a given number of trials, based on the Poisson fluctuations of the signal 


8 


around p 0 {9i,r). 

We define a 5<r detection after correction for Nt trials to be the number 
of counts n which would have a Poisson probability P corresponding to a 
Bonferroni corrected p- value p c given by 

p c = p 0 /N t = P(> n\n hk ). (4) 

Here, po = 2.3 x 10“ ' is the p- value corresponding to 5a, nbk = rR{9i) is 
the number of background counts expected over the burst duration r, and 
P(> n|nbk) denotes the Poisson probability of getting n or more counts when 
the Poisson mean is nbk- We take the value of p, o {0i,r) to be the amount 
of expected signal which would satisfy this criterion 50% of the time. We 
find p 0 ( Oi , r) by estimating the Poisson probability P of finding at least n 
counts to be 50% according to the relation 

P(> n\(n hk + p o (0i,T))) = 0.5. (5) 

By equating the values of p o {0i,T~) found from Equation [5] to p(r,9i,r ) 
in Equation [3] and solving for r, we calculate the maximum distance from 
which a PBH burst could be detected by the high-energy observatory, 


^max(^5 7") — 


I (i-/) 

4:irp o (0i,T ) 



dN{r) 

dE 


A(E, Oi) dE 


( 6 ) 


for various zenith angle bands 9% and burst durations r. Denoting the effec- 
tive field-of-view of the detector for a given zenith band by 


FOV e ff (0j) — 27 t(cos 9 itmin cos0j >max ) sr (7) 

where 0,;, m in and #j )m ax are the minimum and maximum zenith angles in 
band i, the detectable volume is then 

V(r) = £m,r) = (8) 


3.3. Upper Limit Estimation 


If the PBHs are uniformly distributed in the solar neighborhood, the X% 
confidence level upper limit ( ULx ) on the rate density of PBHs bursts (that 
is, the number of bursts occurring locally per unit volume per unit time) 
can be estimated as 


UL X 


m 

VS 


(9) 
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if at the X% confidence level the detector observes zero bursts. Here V is 
the effective detectable volume from which a PBH can be detected, S is the 
total search duration and m is the upper limit on the expected number of 
PBH bursts given that at the X% confidence level zero bursts are observed 
at Earth. Note that for Poisson fluctuations P(0\m) = 1 — X and so m = 
ln(l/(l — X)). If X = 99% and m = In 100 ~ 4.6, the upper limit on the 
PBH burst rate density will be 

4 6 

UL 99 = — . (10) 


4. Results 

4-1. Milagro Limits on the Rate-Density of PBH Bursts 

During the early days of its operation, Milagro had lower angular resolu- 
tion prior to the addition of the outrigger array and used a triggering system 
that did not accept many of the low energy events. Thus, for this search 
we used the last five years of Milagro data: specifically from 03/01/2003 to 
03/01/2008. (Due to various detector-related issues, 7% of the data taken 
during this period was also not used.) Selection cuts were made to increase 
the quality of the data searched. Reconstructed events which have a pre- 
dicted angular reconstruction error greater than 2° were rejected. (This 
corresponds to nfit >20 where nfit is the number of PMTs participating 
in the reconstruction of the shower.) The maximum zenith angle used was 
45° and the best limits were obtained with no gamma-hadron separation cut 
applied, because such a cut also strongly lowered the Milagro photon effec- 
tive area at energies below 1 TeV. Overall, our analysis utilized 1673 days 
(4.58 years) worth of good data, amounting to ~ 93% of the total Milagro 
data collected during the five year period (neglecting the dead time). The 
Milagro search and its optimization presented here are described in further 
detail in m 

We performed a blind search (that is, a search utilizing no external 
triggers) for burst durations ranging from 250 /as to 6 minutes. First we 
created skymaps for overlapping time intervals, each offset by 10% of the pre- 
set burst duration. We then spatially binned the skymap and searched for 
locations with significant excess over background in the Milagro data. The 
optimum bin-size was determined using a Monte-Carlo simulation and varied 
with the pre-set burst duration. For short durations the optimum bin size 
was of order ~1.8° and for long durations it was ~0.8°: for short durations, 
Milagro was signal-limited requiring a larger bin-size to accumulate more 


10 


signal; for long durations, Milagro became background-dominated, requiring 
a more restricted bin-size to reduce background contamination. 

No statistically significant (5a) event was observed over the 4.58 years 
of data. Proceeding on the basis of null detection, we calculated the upper 
limits on the PBH burst rate density following the methodology described 
in Sectional For Milagro, we parameterized the effective area as A(E,6i) = 
10 ai ( lo g£) 2 +o; lo g£+ c ; m 2 w ith the parameters a*, b\ and c« given in Table [l] 
for three zenith angle bands. Figure [2] shows the Milagro effective area 
curves for the selected three zenith angle bands. 


Zenith Angle Band 9i 

di 

bi 

Ci 

0° - 15° (0i) 

-0.4933 

4.7736 

-2.4272 

15° - 30° (0 2 ) 

-0.5037 

5.0102 

-3.4015 

30° - 45° (0 3 ) 

-0.4273 

4.7931 

-4.3030 


Table 1: Milagro effective area parametrization for various zenith angle bands. 


We utilized a Monte Carlo simulation to calculate the /r o (0?., r) values for 
various burst durations. Because the trials in our search were not indepen- 
dent, we took Nt in Equation [4] to be the effective number of independent 
trials calculated using the method described in m- The effective number 
of independent trials ranged from ~0.1% to ~40% of the total number of 
trials depending on the search duration, with the shorter search durations 
having the lower fraction of effective trials. The dependence of the limit 
on the estimated number of independent trials is quite mild ( N?- 018 ) so 

that varying the estimated trials by 3 orders of magnitude produces less 
than 15% change in the limit. The resulting /r o (0j,r) values are listed in 
Table [2] These ^ o (0i,t) values and the Milagro effective area parameteriza- 
tions were then used to derive the maximum detectable distance of a PBH 
burst, r max (0j,r), using Equation [6] assuming an energy range of £i=50 
GeV and i?2=100 TeV and a deadtime of 7%. The derived r max (0j,r) val- 
ues were then used to calculate the effective volume that was probed by the 
Milagro observatory. Using Equation 10, we calculated 99% upper limits 
for various PBH remaining lifetimes and the effective total search period of 
4.58 years. Our results are shown in Tableland in Figure [3] According to 
our results, Milagro is most sensitive to burst durations of about 1 s. For 
shorter durations, the Milagro data is starved for signal photons; for longer 
durations, the background tends to dominate the signal. We note that Mi- 
lagro has a systematic flux uncertainty of ~30% [18] which translates into 
an ~50% uncertainty in our calculated limit (shown as a pink shaded band 
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in Figure [3]). 


Burst Duration r (s) 

Mo ,T 

UL 9g (pc 3 yr *) 

0.001 

11 

3.1 xlO 5 

0.01 

16 

1.2 xlO 5 

0.1 

22 

5.4 xlO 4 

1.0 

35 

3.6 xlO 4 

10.0 

65 

3.8 xlO 4 

100.0 

150 

6.9 xlO 4 


Table 2: Counts (/z 0 (V)) needed over the background for a 5a detection with 99% proba- 
bility and calculated 99% confidence upper limits (ULgg) for various burst durations (r) 
for Milagro. 


4-2. Improved HAWC Sensitivity to the Rate Density of PBH Bursts 

Milagro’s successor HAWC is located at a higher altitude and features a 
better detector design, allowing for superior sensitivity to PBH bursts. In 
this section, we apply our methodology to estimate the HAWC sensitivity to 
PBH bursts. In our calculations, all relevant characteristics of the HAWC 
detector are encoded in the effective area. We calculate the effective area 
using a Monte Carlo simulation which models the interaction of photons 
and cosmic rays in the atmosphere and the response of the detector to the 
extensive air showers generated by these particles. The effective area is then 
defined as the ratio of the number of events that satisfies a given set of 
cuts to the total number of events multiplied by the total throw area of the 
Monte Carlo simulation. In our case the Monte Carlo throw area is a circle 
of 1000m radius. The cuts are comprised of a trigger cut, an angle cut and a 
gamma-hadron separation cut. For the trigger, HAWC will use events with 
nHit, the number of PMTs hit by the air shower, greater than a certain 
value. The angle cut is employed to specify the direction of the photons 
and is a measure of HAWC’s angular resolution. In the simulated events 
we use an angular parameter Del Angle which is the difference between the 
true location of the particle in the sky and the reconstructed sky location 
of the particle. This is a proxy for the angular search bin-size. Because the 
background is predominantly protons, the angle cut was not used to calculate 
the HAWC effective area for protons and the gamma-hadron separation 
cut was used to reduce background events. The standard gamma-hadron 
separation parameter for HAWC is called compactness and is defined as 
nHit/CxPEAO where CxPEA 0 is the number of photoelectrons recorded in 
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Photon Effective Area 



Figure 2: Effective Area for photon detection for HAWC and Milagro as a function of en- 
ergy. The Miiagro effective area curves use nHit >50, DelAngle <1.5 deg and no gamma- 
hadron separation cuts. The HAWC effective area curves use nHit >100, DelAngle <0.8 
deg and nHit/CxPE40 > 7. The HAWC cuts are optimized for the PBH spectrum and 
utilize an nHit cut that is well above the intrinsic threshold. This and the fact that Mila- 
gro used no gamma-hadron cut result in an effective area for Milagro which is larger than 
for HAWC at low energies. However, HAWC has superior sensitivity in the PBH search. 

the strongest hit PMT outside a 40m radius from the reconstructed shower 
core. The shower core represents the location on the ground where the 
original particle would have hit had it not interacted with the atmosphere. 
Figure [2] shows the HAWC effective area A(E,0i) for photons as a function 
of photon energy and zenith angle band using the optimum cuts for a 10 s 
burst search ( nHit >100, DelAngle <0.8 deg, and with nHit/CxPEA 0 > 

7 ). 

In order to estimate the background rate where £ is a measure 

of the spatial resolution of HAWC, we used the ATIC cosmic ray spectrum 
given in Reference m 

ijy / \ — 2.65 

y = 7900 x ( ) particles m” 2 s _1 sr _1 GeV _1 (11) 

dE \ lGev J 

and convolved it with the HAWC proton effective area for a given zenith 
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angle band, 


f E 2 rlAT 

Rb(0i,O= —^ApiE^^dEx 2tt(1.0 - cos(O) x 1.2 (12) 

JEx d E 

In Equation 12, the 27r(1.0 — cos(£)) term represents the patch in the sky 
that corresponds to the spatial resolution of HAWC in steradians. In our 
case £ = DelAngle. The factor 1.2 is a correction factor incorporating the 
other particle species in the cosmic ray background. 


Parameter 

Values 

nHit > 

30, 70, 100, 150, 200, 250 

DelAngle < 

0.1°, 0.3°, 0.8°, 1.0°, 3.0°, 8.0° 

Compactness > 

1, 2, 3, 4, 5, 6, 7, 8, 9, 10 


Table 3: Various parameter cuts used for the simple parameter search. 


Because we seek the sensitivity in the case where there is no prior knowl- 
edge of the burst location, we need to take into account the number of trials 
performed for the search. For example, if the HAWC field-of-view (1.8 sr) 
is divided into spatial bins of radius 0.7°, then there will be approximately 
10 4 spatial bins (trials) per time bin searched. The optimum spatial bin- 
size depends on the search duration, the trigger criteria, and the value of the 
gamnra-hadron separation parameter. The number of time bins is estimated 
by dividing the total search period (estimated as 5 years for HAWC) by the 
burst duration r. Thus the total number of trials depends on the burst du- 
ration r, the optimal spatial bin-size DelAngle , the trigger criterion nHit 
and the value of the compactness parameter. In order to find the optimum 
set of cuts, we performed a simple parameter search and identified the set 
of values which give the best PBH limit according to the method described 
in Section [3l 

For burst durations ranging from 0.001 s to 100 s, we performed cuts 
on all parameter combinations given in Table [3] on the Monte Carlo output 
and calculated corresponding effective areas for photons and protons. Us- 
ing Equations EU and p~2| we then calculated the background rate and the 
background number density n-bk (#*)£) = Ti4(0j,£) (see Equation [4]) which 
depends on the zenith angle band and the spatial resolution. As remarked 
earlier, the effective number of independent trials differs for each parame- 
ter combination. Taking this into account, we have calculated the ^ 0 (#j,r ) 
values corresponding to burst durations ranging from 0.001 s to 100 s for 
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Duration r (s) 

nHit 

DelAngle (deg) 

Compactness 

0.001 

30 

3.0 

3 

0.01 

30 

3.0 

3 

0.1 

70 

1.0 

4 

1.0 

70 

1.0 

5 

10.0 

100 

0.8 

7 

100.0 

100 

0.8 

7 


Table 4: Optimized cuts for various burst durations. 


various zenith angle bands. These /j, o (0i,T) values and the effective area for 
photons are then inserted into Equation [d] (with £l= 50 GeV and £'2=100 
TeV) and the maximum distances r max (0,j,r) at which a PBH burst can be 
detected by the HAWC observatory is calculated for various burst durations 
r assuming negligible dead time. Using these r max (0j,r) values and Equa- 
tion [8j we have determined the effective detectable volume V (r) for each 
burst duration. The PBH limit for each parameter combination was calcu- 
lated using Equation [To] and the set of cuts that gives the best limit for a 
given burst duration was selected. The resulting optimized parameter cut 
values are given in Table [4| The corresponding values of // o (0i,r) for a 5a 
detection are given in Table [5] with the associated background counts and 
number of trials factor. 

The final 99% confidence level upper limits on the PBH burst rate den- 
sity are given in Table [6j assuming zero PBH bursts are observed over the 
5 year HAWC search period. For each burst duration, the maximum de- 
tectable distance for zenith angle band 9\ and the corresponding effective 
detectable volume are also shown in Table [6j We note that HAWC sys- 
tematic uncertainties have not been included in this study. In Figure [3j 
the blue, green, and red thin dashed curves denote the PBH rate density 
upper limits that HAWC will set if zero PBH bursts are detected over a 
one, two and five year search period respectively. Upper limits based on 
earlier null detections from various other observatories are also shown in 
Figure | [201 ED [221 [231 ED]- All limits shown in Figure [3] were obtained 
based on the PBH Standard Emission Model diE]. 
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Figure 3: PBH Burst Rate Density Upper Limits from Milagro and projected for HAWC, 
compared with limits from previous direct search experiments [20i mum mi no]. The 
pink band represents the 50% systematic uncertainty of the Milagro limit. All limits are 
at the 99% Confidence Level (we have rescaled the reported 95% CL HESS limit to 99% 
CL) and were obtained based on the PBH Standard Emission Model. 
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Duration r (s) 

Zenith Angle Band 0i 

Number of Trials 

Bgnd. Counts nbk 

Mo (0i, r) 

0.001 

0° - 26° (6h) 

3.3 x 10 i:i 

0.0637 

10.6 

0.001 

26° - 37° (0 2 ) 

3.3 x 10 i:j 

0.0240 

8.6 

0.001 

37° - 46° (0 3 ) 

3.3 x 10 13 

0.0083 

7.7 

0.001 

46° - 53° (0 4 ) 

3.3 x 10 la 

0.0026 

6.7 

0.01 

0° - 26° (0i ) 

3.3 x 10^ 

0.6372 

17.0 

0.01 

26° - 37° (0 2 ) 

3.3 x 10^ 

0.2397 

13.4 

0.01 

37° - 46° (0 3 ) 

3.3 x 10 12 

0.0832 

10.6 

0.01 

46° - 53° (0 4 ) 

3.3 x 10 12 

0.0256 

8.6 

0.1 

0° - 26° (0 4 ) 

3.0 x 10^ 

0.1355 

11.5 

0.1 

26° - 37° (0 2 ) 

3.0 x 10 rz 

0.0456 

9.6 

0.1 

37° - 46° (0 3 ) 

3.0 x 10 12 

0.0144 

7.7 

0.1 

46° - 53° (0 4 ) 

3.0 x 10 1 ' 2 

0.0036 

6.7 

1.0 

0° - 26° (0 4 ) 

3.0 x 10“ 

1.0481 

18.6 

1.0 

26° - 37° (0 2 ) 

3.0 x 10 11 

0.3422 

14.3 

1.0 

37° - 46° (0 3 ) 

3.0 x 10 11 

0.1055 

10.6 

1.0 

46° - 53° (0 4 ) 

3.0 x 10“ 

0.0251 

8.6 

10.0 

0° - 26° (0 4 ) 

4.6 x 10 iu 

2.4405 

23.2 

10.0 

26° - 37° (0 2 ) 

4.6 x 10 1U 

0.7039 

16.0 

10.0 

37° - 46° (0 3 ) 

4.6 x 10 1U 

0.1912 

11.5 

10.0 

46° - 53° (0 4 ) 

4.6 x 10 iu 

0.0451 

8.6 

100.0 

0° - 26° (0 4 ) 

4.6 x 10 a 

24.4049 

51.3 

100.0 

26° - 37° (0 2 ) 

4.6 x 10 y 

7.0394 

31.6 

100.0 

37° - 46° (0 3 ) 

4.6 x 10 9 

1.9118 

20.8 

100.0 

46° - 53° (0 4 ) 

4.6 x 10 a 

0.4513 

14.2 


Table 5: Counts /j, 0 (di,T ) needed over the background for a 5cr detection with 50% prob- 
ability for various burst durations and zenith angle bands for HAWC. 


Burst Duration r (s) 

r max (pc) 

Effective Volume V(r) (pc 3 ) 

PBH Upper Limit (pc 3 yr x ) 

0.001 

0.033 

0.000016 

56861 

0.01 

0.044 

0.000042 

21976 

0.1 

0.062 

0.000092 

10038 

1.0 

0.078 

0.000172 

5354 

10.0 

0.089 

0.000227 

4059 

100.0 

0.085 

0.000191 

4822 


Table 6: The maximum detectable distance (for zenith band 9\), the detectable effective 
volume and HAWC PBH limit in the event of null detection in 5 years for various remaining 
PBH lifetimes. 
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5. Discussion 


In this work, we report new PBH burst rate density upper limits on par- 
sec scales based on a direct search performed with Milagro data. These new 
Milagro limits probe various burst timescales previously not investigated. 
For 1 s bursts, which were probed by numerous earlier experiments, the 
Milagro limit is now the most constraining. Only the HESS limit at 30 s is 
more constraining than the Milagro limits. Milagro’s successor, the HAWC 
observatory, will be even more sensitive to PBH bursts. As seen in Figure [3j 
a null detection from a 5 year search with the HAWC Observatory will set 
PBH upper limits which are significantly better than the upper limits set by 
any previous PBH burst search including Milagro. Also note that HAWC 
will surpass the current HESS best limit for a 30 s burst in one year. Ac- 
cording to our study, if a PBH explodes within 0.074 parsec (15,000 AU) of 
Earth and within 26 degrees of zenith, HAWC will have a 95% probability 
of detecting it at 5a (as optimized in a 10 s search after trials corrections). 
HAWC would see with 95% probability a PBH burst within 37 degrees of 
zenith if it happens within 0.058 parsec (12,000 AU) of Earth. 

Direct search limits are limits on the local distribution of very low mass 
black holes, regardless of their initial mass or their formation era and mech- 
anism. A 10 9 g black hole has a remaining evaporation lifetime of 1 s and a 
7 x 10 11 g black hole has a remaining evaporation lifetime of 5 years. It can 
be shown that the evaporation rate jT5j of an individual black hole whose 
remaining lifetime is much less than the age of the universe implies that 
the number density function of any population of BHs with present masses 
~ M « 5 x 10 14 g has the form dn/dM oc M 2 around M [23]. From Figure 
3, a null detection from a 5 year HAWC search would correspond to an upper 
limit on the number of local bursts of less than 2.0 X 10 4 pc -3 over 5 years 
and hence to an upper limit on the local density in 7 x 10 11 g or lighter black 
holes of p(< M ) < 5 x 10 _18 Mq pc~ 3 . This applies to very small black holes 
produced in the present universe as well as primordial black holes. This limit 
is well below the total (visible and dark matter) local dynamical mass den- 
sity in the solar neighborhood (the Oort limit) determined from Hipparcos 
satellite measurements, PQ,oort = 0.102(±0.010)Mq pc 3 [25], and the local 
dark matter density in the solar neighborhood, Pq,dm = 0.008(±0.003)Mq 
pc~ 3 [26]. If the present number density function for such small black holes, 
dn/dM oc M 2 , can be extrapolated to black holes of present mass 5 x 10 14 g, 
a null detection from a 5 year HAWC search would correspond to an upper 
limit on the local density in 5 x 10 14 g black holes of p(< M) < 10 _6 Mq 
pc~ 3 , five orders of magnitude less than Po^oort- 
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Direct search limits on the local rate density of PBH bursts are higher 
than the average cosmological PBH burst rate density limit implied by the 
100 MeV extragalactic gamma-ray background constraint on the emission 
of 5 x 10 14 g black holes EU22Q However, PBHs have the gravitational 
properties of cold dark matter and so should be clustered in our Galaxy, 
enhancing the local PBH burst rate density by many orders of magnitude 
over the average cosmological rate [28]. Thus a substantial number of PBHs 
that evaporate producing gamma ray bursts may exist in our Galaxy. PBHs 
clustered in our Galactic halo should also contribute an anisotropic Galactic 
gamnra-ray background, separable from the extragalactic gamma-ray back- 
ground. Wright claims that such a halo background has been detected [5]. 
The direct search limits are also weaker than the average Galactic PBH 
burst rate density limit on kiloparsec scales implied by the Galactic antipro- 
ton background [9] but the antiproton-derived limit depends on the assumed 
PBH distribution, the propagation of antiprotons within the Galaxy, and the 
secondary antiproton component. Direct search limits are independent of 
assumptions concerning the PBH spatial distribution. 

In deriving the new limits, we have assumed the Standard Evapora- 
tion Model in which the Hawking-radiated particles leave the vicinity of the 
black hole without substantially interacting with other Hawking-radiated 
particles. As shown in detail in [3], the energy of the Hawking-radiated par- 
ticles and the time interval between subsequent emissions are such that self- 
interactions between the Hawking-radiated particles or their decay products 
should not form a QED or QCD photosphere around the evaporating black 
hole nor modify the predicted PBH gamma ray burst spectra. However, 
if the PBH is embedded in an ambient high density plasma or strong mag- 
netic field, interactions may arise with the net effect that the predicted PBH 
burst gamma ray spectrum may be enhanced at low energies and decreased 
at high energies. Because of the remaining BH lifetime’s r oc T~ 3 depen- 
dence (see Equation [l]) , as-yet unknown particle modes manifesting at high 
temperatures are unlikely to substantially modify the predicted spectra by 
more than a factor of 2. An exception would be the low temperature Hage- 
dorn Model in which all of the final burst emission is produced at a black 
hole temperature close to the QCD confinement scale, A qcd — 250 — 300 
MeV [7]. The low temperature Hagedorn Model would generate a stronger 


1 We note that all limits on the expected burst rate derived from the extragalactic or 
galactic gamma ray or antiproton backgrounds have assumed an extrapolation of the form 
dn/dM oc M 2 from masses of 5 x 10 14 g down to the very small masses of presently 
bursting PBHs. 
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burst in gamma rays detectable at lower energies but is not consistent with 
the QCD behaviour observed in terrestrial accelerators [3J. 

In conclusion, the HAWC observatory has the ability to directly detect 
emission from nearby PBH bursts. This capability is scientifically very im- 
portant, given the large number of early universe theories that predict PBH 
formation and the uncertainty in the degree to which PBHs may cluster 
locally. A confirmed direct detection of an evaporating PBH would also 
provide unparalleled insight into general relativity and high energy particle 
physics. 
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